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ABSTRACT 

We introduce a prescription for the luminosity from accreting protostars into smoothed parti- 
cle hydrodynamics simulation, and apply the method to simulations of five primordial mini- 
halos generated from cosmological initial conditions. We find that accretion luminosity delays 
fragmentation within the halos, but does not prevent it. In halos that slowly form a low num- 
ber of protostars, the accretion luminosity can reduce the number of fragments that are formed 
before the protostars start ionising their surroundings. However, halos that rapidly form many 
protostars become dominated by dynamical processes, and the effect of accretion luminos- 
ity becomes negligible. Generally the fragmentation found in the halos is highly dependent 
on the initial conditions. Accretion luminosity does not substantially affect the accretion rates 
experienced by the protostars, and is far less important than dynamical interactions, which can 
lead to ejections that effectively terminate the accretion. We find that the accretion rates onto 
the inner regions of the disks (20 AU) around the protostars are highly variable, in contrast to 
the constant or smoothly decreasing accretion rates currently used in models of the pre-main 
sequence evolution of Population III stars. 
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1 INTRODUCTION 

Over the last decade the picture of star formation in the pri- 
mordial universe was that of single massive stars that were the 
sole inhabitants of the first dark matter minihalos formed af- 
ter the big ban g. This was the in f erence made from the numer- 
ical models of I Abel et all feOOOl. |2002[) which was later rein- 
forced by the findings of Yos hida et alj ((2008). In these works 
the star formation process was followed from cosmological ini- 
tial conditions, right up to the densities where the first proto- 
star formed. This required the combination of cutting-edge hy- 
drodynamic simulations with a detailed treatment of the chem- 
istry and thermodynam ics of the gas, as this controls where frag- 



mentation occurs ( e.g. Dalgarno & Lepril 1987t Abel et al 
iGalli & Palli 1 19981; iGlover & Brand! 1200 lb lOmukai et al.l 



1997; 



2005) 



From this lOmukai & P alla (2003) used detailed stellar modelling 
to estimate a final mass of the primordial star and found that above 
an accretion rate of 4 x 10~ 3 M Q yr~' masses in excess of 100 
Mq were produced. 

However the above mentioned simulations were unable to pro- 
ceed beyond the first stage of collapse due to the numerical diffi- 
culty of following the evolution of high density gas, where the nu- 
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merical time step becomes prohibitively small. More recent work, 
which follows the collapse beyond the formation of the first pro- 
tostellar core, has cast doubt on the isolated picture of Population 
III star formation, suggesting that they may have been members of 
binaries, multiples, or even small-N cluster s. Using an idealised 
Bonnor-Ebert model, Machida et al. I J2008h showed that binary 
stars w ere the likely outcome of a rotating minihalo. Idark et all 
d2008h showed that in idealised conditions, gas with a barotropic 
equation of state approximating the behaviour of primordial gas 
could fragment vigorously to form a small cluster. This work took 
advantage of the 'sink' par ticle technique used in present day star 
formation l lBate et al Jl995h to follow the evolution past the forma- 
tion of the first object by replacing high density collapsing gas with 
a non-gaseous particle that can accrete additional bound gas but that 
only interact s with its surrou ndings via gravity. A follow up study 
to this work l lClark et al.l201lb a) found that the fragmentation per- 
sisted when the detailed chemical and thermodynamic evolution of 
the gas was followed. 

Although these simulations make use of idealised initial con- 
ditions, other work has addressed this issue with cosmo l ogical ini- 
tial conditions and full chemical networks. iTurk et al] d2009l) us- 
ing the AMR (adaptive mesh refinement) method showed that a 
wide binary with a separation of 800 AU formed i n one out of 
the five minihalos that they simulated. IStacv et alj ( |2010|) used 
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smoothed particle hydrodynamics (SPH) combined with the sink 
particle technique to show that there is further fragmentation after 
the first object for ms in a minihalo an d that a small multiple sys- 
tem can be formed. Idark et alj ( 1201 ll b) found that the individual 
disks around Population III stars are prone to fragmentation and are 
likely to fragment into h igher-order multip les. This result was con- 
firmed more recentl y bvlGreif et al.lfcOl ll) using the novel moving 
mesh AREPO code 1 Sprin gell20icl) . 

The case for fragmentation therefore seems robust, as it has 
been found by multiple authors using complementary numerical 
techniques. If this result withstands further investigation, then it 
will have important implications for our understanding of cosmol- 
ogy and the early universe. For instance, the first stars are thought 
to be an important source of ionizing photons in the early uni- 
verse, and hence an i mportant contributor to the reionisation of the 
intergalactic medium dKitavama et al". 2004; So kasian et alj|2004l ; 
IWhalen et al.l2004l ; lAlvarez et al]2006t I Johnson et al.l2007t) . If the 
mass available for accretion is split into multiple stars there may be 
fewer truly massive stars, which would reduce the number of ioniz- 
ing photons produced. Likewise, the eventual fate of a Population 
III star, and consequent ly its enrichment of the surrounding gas, is 
determined by its m ass dFrver & KalogerakOOltlHeger et al.l2003l : 
lYoshida et alj|2004h . Additionally, while the traditional picture of 
extremely massive single stars meant there would be few observa- 
tional signatures that would be observable today, the new picture of 
a few multiple stars l eads to some new observational possibilities. 
In the simulations of iGreif et ail J201 lh there are many dynamical 
interactions between the protostars that lead to the ejection of some 
of the low-mass protostars. There is therefore the exciting possibil- 
ity that these stars could still be shining today, providing a direct 
insight into the physical conditions in the high-redshift universe. 
Moreover the p ossibility of tight binary systems resulting from disk 
fragmentation dClark et aill201 ll b) would allow primordial X-ray 
binaries, cataclysmic variables, or even gamma ray bursts. 

Given the potential implications of these results, it is vital 
to consider whether there are any mechanisms that could sup- 
press fragmentation within primordial minihalos. Ionisation from 
the protost ar is the most likely mechanism to stop further frag- 
mentation. Om ukai & Inutsukal |2002j) showed that above a criti- 
cal ionising flux an HII region from a primordial star can unbind 
the surrounding gas, while ionising fl uxes below th i s valu e may 
nevertheless suppress fragmentation. Tan & McKee (2004) use a 
semi-analytic model to describe the evolution of feedback from a 
steadily growing protostar and find that once the star is older than 
its Kelvin-Helmholtz time, and is contracting towards the main se- 
quence, there is a rapid increase in the amount of ionising radiation 
it emits. One could speculate that this transition may therefore mark 
the point at which fragmentation is halted. In simulations of local 
star formation, ionisation has bee n found to be unable to prevent 
fragmentation jPeters et al.ll2010l a,b,c). However, it is likely that 
the effect of ionising radiation was more significant in the early uni- 
verse, since it would have been accompanied by the dissociation of 
H2 by Lyman- Werner photons, thereby removing the primary gas 
coolant iomukai & Nishill999l ; lGlover & BrancfeOQll) . 

The early period of protostellar growth before the ionising ra- 
diation becomes important therefore represents the most favourable 
window of opportunity in which fragmentation can occur. One of 
the few processes that would act against fragmentation in this epoch 
is accretion luminosity feedback from the forming pro tostars. We 
first in troduced this in the simulations reported on in IClark et all 
( l201ll . b) in order to study its effects on the stability of Pop. Ill 
accretion disks. We found that accretion luminosity does indeed 



change the disk evolution, but that it cannot ultimately prevent the 
disk from fragmenting. The feedback is able to support the inner 20 
AU of the disk, which was previously unstable, against fragmenta- 
tion, but the outer regions still fragment, albeit after a longer time 
period. A more detailed desc ription of how fra gmentation in the 
disk proceeds can be found in IClarketal.l l l20TTl b). 

In this work we seek to address a different question. The opti- 
mal time for fragmentation within a minihalo is the first few thou- 
sand years before the first protostar approaches the main sequence, 
at which point ionisation feedback will become significant and will 
act against further fragmentation. In this work we aim to capture 
the full evolution of the halo during this regime to determine how 
much fragmentation can occur in this time interval, and to what 
extent radiation affects the fragmentation. 



2 THE METHOD 

We perform the calculation s for this paper using the SPH code 
GADGET 2 dSpringeJl2005h . We have substantially modified this 
code to include a fully time-dependent c hemical network, d etails of 
which can be found in the appendix of Clark et alj (1201 lL a). Our 
treatment includes: H2 cooling dGlover & Abel 200 8]); optically 
thick H2 cooling using the Sobolov approximation ( Yoshida et al. 
2006); collisionally induced H2 emission {Ri pamonti & Abe] 
2004); ionisation and recombination dGlover & Mac Lowl2007 ): as 
well as heating and cooling from changes in the chemical makeup 
of the gas, and heatin g and cooling from shocks, compression and 
expansion of the gas. lTurk et al.l d201 ll) showed that the choice of 
H2 three-body formation rate coefficient influences the resulting 
dynamics of the gas within t he halo . In this work we use the three- 
body H2 formation rate of iGloverl (|2008) which is intermediate 
within the range of the published rates and is based on applying 
the principle of detailed balance to a relatively recent calculation of 
the collisional dissociation of H2. 

We include heating from the accretion luminosity as an ad- 
ditional heating term when solving the ordinary differential equa- 
tions that govern the chemical and thermal evolution of the gas. 
While the protostar will also have an interior luminosity, for the 
majority of the early protostellar evolution t his is an order of mag- 
nitud e lower than the accretion luminosity dHosokawa & Omukail 
2009), and so we focus here on only the accretion luminosity. The 
accretion luminosity will also transfer some momentum to the gas. 
However, the resulting outward force is many orders of magnitude 
smaller than that of the gravitational force on the gas, and so it is 
safe to neglect it. 

The accretion luminosity is calculated from the standard equa- 
tion, 



GM t M 



(1) 



where M is the accretion rate of the protostar and R„ is the stellar 
radius. We make the assumption that the gas is optically thin to 
the emitted radiation which ensures that we are overestimating the 
effects of the accretion luminosity to obtain an upper limit on the 
feedback effects. The heating rate per unit mass for the gas will 
then become 
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where p g is the gas density, Kp is the Plank mean opacity and 
r is the distance of the gas from the source. We calculate the 
Planck mean opacity of the gas by interpolating from the tables 
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of iMaver & Duschll (2005) which include contributions from deu- 
terium and lithium in the gas in addition to hydrogen and helium. 

To accurately calculate the accretion luminosity, both the ac- 
cretion rate and the stella r radius need to b e known. We achieve this 
by using sink particles dBate et aT] 1 19951) to model the protostars 
and record their growth in mass throug hout the simu l ation. These 
were first implemented into Gadget bv ljappsen et al.l l| 2005i) . Sink 
particles are non-gaseous particles that replace extremely dense gas 
if it is found to be both gravitationally bound and collapsing. This 
allows us to evolve of the simulation beyond the point at which 
the first object form s. For a recent d i scussi on of the pro's and cons 
of sink particles see lFederrath et all j2010t) . We form sink particles 
above densities of 10 15 cm~ 3 , after which there are no more chem- 
ical heating terms that can prevent the gas collapsing to form a pro- 
tostar. The sinks have accretion radii of 20 AU and consequently the 
inner accretion disk is not resolved (alt hough an y disk outside this 
distance is resolved). Iciark et al.l d201lL b) and iGreif et ai] boilh 
consider fragmentation within this regime and the effect of accre- 
tion luminosity upon it. In this work, we use larger sink radii to 
allow us to study the cluster as a whole. 

The accretion rate of the sinks can be found from simply look- 
ing at how their mass grows in time. However, as SPH is a particle- 
based method, accretion occurs in discrete units and can be noisy. 
In order to account for this, we calculate the accretion rate by tak- 
ing a smoothed average over the last 100 years, updated at 10 year 
intervals. For the typical accretion rates of the sinks in our simula- 
tion this equates to between 0.1 and 1 M© of accreted material. This 
is equivalent to the accretion rate being smoothed over 10 3 - 10 4 
particles (and updated only after at least 100 new particles have 
been added) and therefore we can be sure that any variability in the 
accretion rates is not due to particle noise. 

In reality, material will flow on to the protostar through the 
inner disk, which will delay it from reaching the protostar. Inwards 
mass transport typically takes place over the viscous timescale, 
which for a thin disk i s typically much larger than the dynamical 
timescale lPring ldll98lh . Howe ver in the primordial protostellar 
disk study of IClark et al.l (20U b), the disk had a scale height of 
approximately 5 AU which, given that we are considering an inner 
disk of 20 AU, means that the thin disk approximation is not valid. 
Furthermore, the disk was self-gravitating and was rapidly trans- 
ferring its angular momentum through gravitational torques which 
lead to high accretion rates. Given these findings, as we do not re- 
solve this region, we simply update the accretion rate immediately 
after it is calculated. However, the procedure of averaging the ac- 
cretion rates which we adopt for numerical reasons will to some 
extent mimic the effect of accreted material being buffered by the 
inner disk. 

Accurately finding the stellar radii is a complex problem 
that would require the implementation of detailed stellar evolu- 
tion models within our hydrodynamic simulation which is beyond 
current computa t ional resources. Instead we used the models of 
lOmukai & Palial J2003h to derive a simple power law approxima- 
tion of the stellar radius. In the early stages of the protostellar evo- 
lution the cooling tim e of the interior is ve ry long and the protostar 
evolves adiabatically. Stahler et al. I Jl986l) showed that the stellar 
radius during this phase grows according to 
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where is the stellar radius in R0, M* is the current stellar mass 
in M© and M is the accretion rate in M@ yr . After some time, 
the internal heat becomes sufficient to drive an outward luminosity 
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Figure 1. The st ellar radius as a functio n of mass found in the stellar evo- 
lution models of lOmukai & Pallal |2003|) for three different accretion rates. 
The dotted black lines show the stellar radius given by our semi-analytic 
model for these accretion rates. 



wave, which results in the rapid swelling of the stellar radius. Once 
the luminosity wave reaches the stellar surface, the interior achieves 
a relaxed state and undergoes Kelvin-Helmholtz contraction until 
the radius shrinks to its main sequence value. 

Fi gure[[]shows tne variatio n of stellar radius with mass calcu- 
lated by Omuk ai & Pallal ( 120031) . Using the model with an accretion 
rate of 10~ 3 Mq yr~'as our fiducial case, we found that the stellar 
radius for this accretion rate could be described by, 
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where p\ = 5 M© is the transition between the adiabatic phase and 
the luminosity wave, and p2 = 7 M© is the transition between the 
luminosity wave phase and the Kelvin-Helmholtz stage. The con- 
stants Ai and A2 are determined at these transition zones to give a 
continuous function. At R m the radius has shrunk to its main se- 
quence value and the feedback will be totally dominated by ionis- 
ing radiation. The main sequence radius in Rq for these stars from 
lOmukai & PaiiaU2003h is 

R ms = 0.28M° 61 R Q (5) 

To generalize this prescription to the case where M 7^ 10~ 3 
Mq yr~ 1 , one must account for the fact that the transition points 
Pi and P2 between the phases scale with the accretion rate as 



-5M' 



0.27 
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(6) 



The constants A 1 and A2 must be correspondingly adjusted to main- 
tain a smooth function. The resulting model is over-plotted onto 
Figure[T]and can be seen to qualitatively capture the change in stel- 
lar radiu£_widijrmss ; _An important caveat here is that the models 
of lOmukai~& Palla ( 2003) are calculated using a constant accretion 
rate, which we shall later to show to be a poor assumption. How- 
ever, the current model is a first step that couples the accretion rates, 
masses and stellar radii in real time with both the dynamics and the 
effects of the protostar's own radiation. 



4 Smith et al. 




10 



100 



1000 



10000 



Column Density [g cm ] 



Figure 2. Column density projection of the centre of Halo 1 just after the first burst of fragmentation simulated at our standard resolution and at ten times 
higher resolution. Feedback is present in both cases. The same number of fragments form in each case and in similar places. 



3 INITIAL CONDITIONS 

In order to accurately capture the true properties of the first mini- 
halos within w hich primordial sta rs form, we use the co smological 
simula tions of iGreif et al.l d201 ll) as initial conditions. iGreif et all 
presented simulations of primordial minihalos that strongly 
fragmented. T hese simulatio ns made use of the novel moving mesh 
code AREPO ( Springell uOlOl) to fully resolve the formation of five 
minihalos from truly cosmological simulations. Cells were refined 
during the evolution to ensure that the Jeans length was always re- 
solved by at least 128 mesh points. The refinement was deactivated 
at densities greater than tin = 10 9 cm~ 3 by which point the mass of 
each element was around 10~ 4 Mq . All of the halos modelled by 
Grei f et al. I J20T 1) form multiple protostars with a range of masses. 
Some protostars are even ejected and stop accreting entirely. These 
minihalos, therefore, represent the perfect sample in which to test 
the effect of accretion luminosity upon fragmentation and to com- 
pare the magnitude of any effects to those of cosmic variance and 
dynamical interactions. 

For th is work we cut out the central two parsecs of the 
IGreif et all simulations and continue their evolution using our mod- 
ified version of Gadget 2 with feedback as discussed in the previous 
section. Each mesh point in AREPO is converted to an SPH particle 
with the same properties as the original element. As the chemical 
network that is implemented in AREPO is based on the network we 
use here dClark et alj|201 3, b) the chemical abundances can also 
be transferred from the original cosmological simulation. All the 
fragmentation and accretion takes place in the central region of the 
halos where the SPH particle masses a re 10~ 4 M fl which gi ves us 
a mass resolution of at least 10~ 2 Mq teate&Burkertlll997l) . 

Table [T] shows the initial conditions of the five halos at the 
point where our simulations commenced. Each halo is simulated 
twice, once with and once without feedback as a reference case. 
The halos typically contain a few thousand solar masses, with the 
largest having 3000 Mq . The mean densities are over 10 7 cm~ 3 
and the mean temperature are around 300 — 500 K. The gas is pri- 
marily atomic with most of the molecular hydrogen being con- 
tained within a disk-like region at the centre of the halos. A more 



Table 1. The initial state of the inner 2 pc of the live minihalos. From each 
minihalo two simulati ons were run, one with feedback, and a control run 
with no feedback. See lGreif et at] fepl ll) for a more detailed description of 
the halos. 



Name M [ Mq ] n [ cm" 3 ] T [K] 



mhl 


1810 


4.62 x 10 7 


409 


mh2 


1240 


8.09 x 10 7 


329 


mh3 


1030 


6.31 x 10 7 


292 


mh4 


2000 


7.92 x 10 7 


458 


mh5 


3340 


8.60 x 10 7 


494 



detailed discu ssion of the physic al properties of the considered ha- 
los is given in lGreif et al.l d201ll) . 

As a resolution study we also increas ed the resolution of mhl 
by a factor of ten using particle splitting dKitsionas & Whitworthl 
|2002|) . Figure [2] shows a column density projection of the centre of 
halo 1 with feedback in our standard run and in the increased res- 
olution run. In the higher resolution run smaller scale structure is 
resolved, but the same number of sinks are formed in both cases. 
The increased resolution decreases the numerical viscosity of the 
simulation which allows the disk to form earlier in the simulation 
as it drains more slowly due to the reduced angular momentum 
transport. Consequently the high resolution image of Figure [2] is 
shown at an earlier snapshot in the simulation than the low reso- 
lution case. However the fragmentation of the disk is qualitatively 
the same once it forms in both cases. The relatively large accre- 
tion radius of our sink particles (~ 20AU) means that all structure 
smaller than this scale is swallowed by the sink and this limits the 
differences between the high and low resolution runs. 



4 RESULTS 

4.1 Fragmentation 

Figure [3] shows a column density projection of the central regions 
of halos 1 and 4 which are good examples of the fragmentation 
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Figure 3. Column density projection of fragmentation seen in the central 1300AU in the first few thousand years for halos mhl (top) and mh4 (bottom). Sink 
particles are denoted by yellow squares. Despite the presence of accretion luminosity heating there is still fragmentation. 



seen in all of the halos. In each case, a disk-like structure is formed 
due to the inability of the halo to transfer angular momentum out- 
wards quickly enough during collapse. While the central region is 
disk-like, it is more extended and irregular than a classical disk. 
In every case this region fragments. Generally several fragments 
form almost simultaneously as the conditions for fragmentation are 
reached at multiple locations within the disk. Halos 1 and 5 frag- 
ment vigorously, whereas halos 2 and 3 fragment more slowly. Halo 
4 is the case that is most affected by accretion luminosity and has a 
rate of fragmentation intermediate between the other cases. 

Let us now consider whether the fragmentation seen here is 
resolved. While the central regions of our halos are disk-like, they 
are strongly self gravitating and have a large vertical extent. As an 
example let us consider the central regions of halo 4 which bears 
the greatest resemblance to a classical disk out of the set of simula- 
tions. Taking a density cut of material greater than 10 ll) cm~ 3 the 
disk-like region had a radius which varies between 200 to 300 AU, 



and has a vertical extent of 150 AU measured from the mid plane. 
The thin disk approximation requires that the vertical extent of the 
disk, H, is much smaller than the radial extent of the disk, R, and 
is onl y approp r iate fo r disks where H/R ~ 0.1 or lower dLodatol 
2008). iNelsonl J200fj|) propose that to avoid artificial fragmenta- 
tion in disks the scale height must be resolved by 4 SPH particle 
smoothing lengths per scale height. In the inner 100 AU of the 
disk in halo 4, the average particle smoothing length is 14.3 AU. 
However, it is unclear how best to define a scale height for this ir- 
regular puffy disk. Given that even our best-case scenario cannot 
be considered a classical disk, and that as the simulations proceed 
fragmentation increasingly takes place in large spiral arm features 
or filaments (as seen in the later panels of Figure [3), a better reso- 
lution criteria to use is the Jeans mass. At our sink creation density 
of 10 15 cm~ 3 , the Jeans mass is 2 x 10~ 2 M e , at least two times 
more massive than our resolution limit. Further to this, as shown 
in Figure [2] even when our resolution was increased by a factor of 
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Figure 4. The time of formation of sink particles and their growth in mass as 
measured in halo 4 without feedback. Each line represents the growth of an 
individual sink. The lower panel shows the standard Gadget simulation used 
here and the upper panel shows a comparison simulation using AREPO with 
the same sink properties. The method of simulation makes no significant 
difference to the number of fragments formed and their growth in mass in 
the overlapping time period. 



ten, the same number of fragments were formed. In truth, the great- 
est limitation in our treatment of the disks is that the large size of 
our sinks means that we do not fully capture the behaviour of the 
inner regions of the disk. As such, there would probably be more 
fragmentation in the disks in reality than we find in our simulation. 
For a fuller discussion of the physics of disk fr agmentation around 
population III protostars see dClark et al1l201 ll b) where this issue 
is considered in great detail. 

The fragme ntation here differs slightly from that shown in 
iGreif et al.l 0oTTh . as in that case the sink radius was similar to 
the stellar radius , whereas here it is 20 AU. The simulations of 
IGreif et al.lj2oT 1) used the AREPO method which has less artificial 
viscosity than Gadget. In order to test whether the results here were 
reproducible with AREPO, we re-simulated halo 4 (which showed 
the maximum difference between the feedback and no feedback 
cases) using AREPO with no feedback and the same sink sizes 
used in this study. Figure [4] shows the growth of fragments formed 
in both simulations at the overlapping times. There is a similar in- 
terval in both cases between the first sink forming and the first burst 
of fragmentation. In both cases the same number of sinks form, al- 
though there are slight differences in the mass growth rates due to 
the different N-body dynamics which occur in each simulation. As 
our results are reproduced by two highly complementary numerical 
schemes, we are confident that we capture the true physical evolu- 
tion and are not strongly influenced by numerical artefacts. 

The fact that larg er sinks are used here compared to the orig- 
inal [Greifetal] d201 ll) simulations mean that we are not resolving 
tight binaries and missing some young low-mass objects formed 
within this radius that may have been ejected. Therefore our 20 
AU sinks are a conservative estimate of the level of fragmentation. 



However, at this radius we are avoiding many of the uncertain- 
ties associated with protostellar mergers. As our young protostars 
would actually be puffy extended objects with radii about 100 Rq 
dStahler et al.ll986l) . there will be strong tidal forces evoked during 
close interactions, leading to the possibility that fragments formed 
close to each other will merge when they interact. It is still un- 
clear how best to treat this possibility. By not forming low- mass ob- 
jects in close proximity to existing sinks, encounters that are close 
enough for the stellar radii to touch occur rarely compared to the 
original simulations (typically between 0-2 times in each halo) and 
we generally avoid this issue. However, despite these small differ- 
ences, qual i tative ly the evolution of the halos is similar to that in 
IGreif et all l l201ll) . with the main difference being that we follow 
the evolution for ten thousand years compared to the original thou- 
sand. 

4.2 The effect of accretion luminosity 

Figure[5]shows the combined mass function of all the sinks formed 
in minihalos 1-5, one and two thousand years after the first sink 
formed. At 2000 yr in the non-feedback case ionisation effects are 
becoming important within Halo 5. However, for the sake of the 
mass function only, we run Halo 5 until this point despite the lack 
of ionisation in our model. This is due to the difficulty in achieving 
a statistically significant number of sinks for the mass function. At 
these early stages the sinks represent protostars rather than finished 
stars, and so these masses will not be those of the final population 
III stars. Nonetheless it can already be seen that the resulting mass 
function will contain a range of masses rather than just being one 
characteristic mass. The mass functions show no systemic varia- 
tion between the case with feedback and the reference case, and 
both cases contain a similar total amount of mass in stars at each 
time. Hence the feedback has not significantly altered the fragmen- 
tation and mass growth when considering the five minihalos com- 
bined. This suggests th at the results of p revious studies which ne- 
glected this effect (e.g. IStacv et aT1l2010h will still be broadly cor- 
rect. The mass functions ap pear to be flatter than the IM F's seen in 
the present day universe (Kroupa 20021. IChabrienl2003b . although 
as yet we have only of order ~ 50 sinks, so this remains statistically 
uncertain. 

Tables[2]and[3]show the number of fragments formed in each 
halo when the mass of the most ma ssive protostar first re aches 
10 or 15 solar masses, respectively. iTan & McKeej d2004l) find 
that ionising feedback does not become effective until the star is 
older than its Kelvin-Helmholtz time and is contracting towards 
the main sequence. For their fiducial model this equates to a 
mass of around 30 Mq for a rotating protostar. However the ac- 
cretion rate for the most massive object is typically only a few 
10~ 3 Mq yr _1 when the pr otostar has reached 10 Mq in our miniha- 
los, whereas in the fiducial lTan & McKed d2004l) models the accre- 
tion rate is 10~ 2 Mq yr~' for a 10 Mq protostar. Since the Kelvin- 
Helmholtz contraction stage commences earlier with a lower accre- 
tion rate, as shown in FigureQ] we estimate that ionisation feedback 
will become important for our minihalos when the most massive 
star is between 10 — 15 Mq . H2 photodissociation will also be- 
come important at this time. Beyond this point, the assumptions 
that we make for the luminosity heating model break down, so we 
chose to terminate the simulations here. 

Perhaps the most striking feature of our calculations is the 
amount of variability between the outcomes of the different halos, 
as shown in Tables [2] and [3] Halo 5 forms a star greater than 10 
Mq after only 600 yr without feedback, whereas Halo 3 takes over 
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Figure 5. The combined mass function of the sink particles formed in the five minihalos one and two thousand years after the first sink particle forms. The 
solid red lines shows the mass function of the halos with feedback, and the dotted black line the halos in the reference case without. In both cases the mass 
function is flatter than the slope of the Salpeter IMF shown by a solid black line. 



Table 2. The number of stars and time when the most massive star in 
the minihalo had a mass of 10 M . . The plus sign next to the number of 
fragments for the halo 3 reference run indicates that no star reached 10 
Mq before the simulation was ended, meaning that the number of fragments 
in this case is a lower limit. The effects of ionising radiation are expected to 
become important once at least one star has reached a mass of 10 — 15 Mq or 
greater, and are likely to suppress further fragmentation. There is signifi- 
cant inter-halo variation both in the number of fragments and the duration 
over which accretion luminosity is the dominant feedback mechanism. For 
equivalent halos, the one which forms the massive star most quickly has the 
least fragmentation. 



Halo 


Ref. 




Feedback 






No. of stars 


Time [yr] 


No. of stars 


Time [yr] 


1 


10 


1,520 


10 


2,520 


2 


10 


7,640 


7 


4,490 


3 


5+ 


9,430 


5 


5,140 


4 


17 


7,320 


5 


1,010 


5 


7 


604 


18 


1,440 



Table 3. The number of stars and time when the most massive star in the 
minihalo had a mass of 15 Mq . A plus sign next to the number of frag- 
ments denotes where a mass of 15 M© was not achieved before the end of 
the simulation and as such the number of fragments is a lower limit. The 
results are similar to Table [2] but the times are longer and there is more 
fragmentation. 



Halo 


Ref. 




Feedback 






No. of stars 


Time [yr] 


No. of stars 


Time [yr] 


1 


11 


2,910 


16 


6,040 


2 


13 


15,020 


8+ 


10,000 


3 


5+ 


9,430 


6 


11,270 


4 


20+ 


22,360 


6 


3,700 


5 


17 


1,060 


23 


3,900 



10,000 yr to do so. In the time to taken to form a 10 Mq protostar, 
the number of fragments varied between 5 and 17 between the ha- 
los. Therefore inter-halo variability is at least as important an effect 
as accretion luminosity feedback. The variability of the halos can 
be traced back to their chemical evolution during their collapse. 
iGreif et al.1 fcoilh found that in two of the halos considered here 
(halos 2 and 3) there was significant HD cooli ng which allowed 
the gas to cool to temp eratures as low as 100K ( Ripamont ij |2007l : 
iMcGreer & Brvanll2008l) . When the gas was reheated by compres- 
sion in the final stages of the collapse this smoothed out some of the 
small scale structure, resulting in less fragmentation. In the remain- 
ing three halos, HD cooling was not activated and so temperatures 
only as low as 200K were obtained via H2 cooling. In this case, the 
subsequent reheating was less violent and more small scale struc- 
ture was retained. A similar reduction in fragmentation has been 
seen in simulatio ns of Pop III. 2 star formation due to reheating 
dClarketalfeOlll ai. 

The first panel of Figure [6] shows the evolution of fragmenta- 
tion within the halos. Halos 1 and 5 fragment rapidly, meaning they 
quickly become dominated by chaotic dynamical interactions and 
the reference and feedback cases are no longer equivalent. Dynam- 
ical interactions are therefore as important as accretion luminos- 
ity effects in halos that fragment rapidly. For example, in the run 
of Halo 1 with feedback, there was a dynamical interaction which 
ejected the most massive star before it could reach 10-15 Mq , and 
consequently there was more time for fragmentation until one of 
the originally lower mass objects reached this mass. Such are the 
numbers of sinks formed within Halo 5 that chaotic N-body inter- 
actions cause the feedback and reference cases to swiftly diverge. 
Consequently no clear conclusions can be made about the effect of 
feedback in Halo 5 and its evolution is not shown in Figure|6] 

Halos 2 and 3 are more straightforward as these halos frag- 
ment and accrete material less vigorously. As fewer fragments are 
formed, there is less competition to accrete the gas, which allows 
the first fragments to grow in mass and substantially heat their sur- 
roundings. This delays when the fragmentation occurs in the feed- 
back case compared to the reference case. In Figure [6] fragmenta- 
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Figure 6. The combined properties of the sinks formed in each halo plotted against the time after the first fragment forms. The solid red line shows the case 
with feedback and the black dotted line the case without. The vertical solid red and black lines in panel one show the point at which a star reaches 15 Mq and 
the dotted vertical lines in panel two show when 10 M ., has been reached. 



tion in the feedback case generally lags behind that in the reference 
case, and in some cases the delay can be as great as a thousand 
years. This was also true in Halo 1 until the chaotic dynamics made 
the runs diverge after a thousand years or so. The delaying of frag- 
mentation seems to be the major consequence of accretion luminos- 
ity feedback. Th is was also the con clusion reached in the protostel- 
lar disk study of IClark et al.l d20 1 ll b). Although feedback does not 
prevent fragmentation, the delay means that there are fewer frag- 
ments when ionising feedback becomes effective, so in total, ac- 
cretion feedback has reduced the number of protostars formed. 

Halo 4 is the case in which the maximal effect of the feed- 
back was seen. As in Halos 1 and 5, HD cooling was not activated 
and 17 protostars were formed in the reference case. However, with 
feedback the number of fragments formed before one of the proto- 
stars reached 10 was reduced from 17 fragments to 7. In Halo 
4 there was enough delay before the second bout of fragmentation, 



after the disk first became gravitationally unstable, that the sinks 
were massive enough to produce a large amount of luminosity. This 
effect was enhanced by the geometry of the resulting system (as 
seen in Figure [3} with all the sink particles remaining within the 
central disk-like region where they could heat the dense gas. 

Panel two of Figure [6] shows that the total mass that goes into 
sink particles shows no clear correlation with feedback. Panel three 
of Figure [6] shows that the total accretion rates have considerable 
temporal variation. This is firstly due to the fact that every time a 
new sink is formed it rapidly swallows up the gas that is bound to 
it, adding a large contribution to the total accretion rate. Once this 
gas has been accreted there is less available for the other sinks and 
the accretion rate c an fall, a process w hich we term fragmentation- 
induced starvation dPeters et al.l2010l c). Another contributing fac- 
tor to the total accretion rate variability is the effect of sink interac- 
tions, which we will discuss more fully in section 1431 In both the 
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Figure 7. The gas heating rates in Halo 1 at a typical snapshot. Red shows 
compressive heating, green shows heating due to H2 formation, and blue 
shows the heating due to accretion. The heating rate from accretion is of the 
same order as that from H2 formation, and about two orders of magnitude 
lower than that from compression during the collapse. 



1CT 5 




n [cm 3 ] 

Figure 8. The gas cooling rates in Halo 1 at the same time as Figure[7j Blue 
shows H2 line cooling, green shows cooling due to collisionaly induced 
emission, red shows cooling from re-expansion of previously compressed 
gas, and pink shows the cooling due to H2 dissociation. 



feedback and reference cases the total accretion rates are broadly 
similar and have a value in the vicinity of 10~ 2 Mq yr~' initially, 
decreasing thereafter. The fourth panel of Figure |6]shows the total 
luminosity of the sinks as a function of time. The total luminosity 
output of the stars has a value of over 10 4 Lq for the vast majority 
of the time. Fragmentation was suppressed most effectively in Halo 
4 and the luminosity shown in Figure|6]shows that this was indeed 
the case in which the feedback was most significant. 

To understand why the accretion luminosity was unable to pre- 
vent fragmentation entirely, we need to consider the chemistry and 
thermodynamics of the gas more fully. Figure UJ shows the heating 
rates experienced by the gas from the various heat sources in Halo 
1 in which feedback was largely ineffective. The heating rate from 
accretion luminosity is two orders of magnitude less than that from 
compression, and about equal to that from H2 formation. We cal- 
culate the heating rate from compressional heating using the below 
formula that can be derived from energy conservation 

dz 

= -e Y V-v (7) 

at 

where £ is the thermal energy per unit volume, yis the adiabatic in- 
dex of the gas and v is the velocity. The heating from compression 
and H2 formation was already being balanced by cooling from H2 
line emission, and re-expansion of the gas, as shown in Figure [8] 
Therefore, the addition of accretion luminosity feedback represents 
only a small change in the thermodynamic equilibrium of the halo. 
However, an important qualification is that luminosity feedback is 
an effect that varies with position, i.e. it is most effective close to 
the sinks where the gas is densest and fragmentation occurs. Addi- 
tionally, the extra heating increases the collisional dissociation rate 
of H2, making it harder for the dense gas to cool. Consequently, the 
accretion luminosity heating is more dynamically significant than a 
first glance at Figure UJ would suggest, which is why it was able to 
delay fragmentation in the minihalos. 



4.3 Accretion and Dynamics 

The previous section demonstrates that accretion luminosity can 
affect the fragmentation seen in the minihalos. However it has also 
shown that these effects can be completely masked by the dynamics 
of the gas. A similar conclusion is derived from massive star forma- 
tion calculations in the present d ay, where dynamical effects dom- 
inate over radiative feedback dKrumholz et al.l [20091 ; IPeters et"al] 
|2010L a,b). To explore this more fully, in this section we contrast 
dynamical effects with those of accretion luminosity upon the in- 
dividual protostars, and consider how this will affect our feedback 
model. 

Figure [9] shows the evolution of the sinks formed in halos 1- 
4. The first sink forms at the centre of the disk and quickly grows 
in ma ss with a smoothly decreasing accretion rate. lYoshida et al.l 
( 2006) showed that the expected accretion rates should be as high 
as 10 _1 M s yr~' after the first half solar mass has collapsed, and 
then smoothly decrease to a value of of 10~ 3 Mq yr _1 when 100 
Mq has collapsed . Our accretion rate for the initial sink agrees 
with lYoshida et all d2006h until fragmentation sets in. At this point 
the accretion rate may briefly rise as the portion of the disk between 
the original sink and the new sinks is strongly torqued, resulting in 
a large outward transfer of angular momentum and inflow of gas 
onto the central sink. After this short transient, the accretion rate 
decreases as the mass available for accretion is now shared between 
multiple sinks. 

Once multiple sinks are formed, the accretion rates of the sink 
particles become highly variable. Due to the high densities charac- 
teristic of primordial star formation, the Jeans length is extremely 
short, and therefore fragments are formed close to each other, lead- 
ing to interactions on timescales comparable to the local free-fall 
time. Figure |9]also shows the paths of the stars in the central 4000 
AU during the period studied here. The sinks that remain in the halo 
centre orbit each other, leading to a periodically varying accretion 
rate as they move around in the gas bound to them. Mor eover, ejec- 
tions a re common and occur in every halo, as shown bv lGreif et al.l 

1). Surprisingly, even the originally central star can be ejected 
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if it has a close three-body interaction as shown in Halo 1 (this only 
occurred in the case with feedback, which is why it was hard to 
compare Halo 1 with and without feedback). Indeed, while we find 
that feedback from accretion luminosity has no significant effect on 
the accretion rates of the protostars, dynamical interactions are ex- 
tremely effective at halting an ejected proto star's accretion entirely 
dReipurfh & ClarkekOOlHBate et alj2002h . 

The accretion rates for the sink particles are variable, yet the 
stellar radius model used to calculate the accretion luminosity was 
developed from simulations with a constant accretion rate. Figure 
1101 shows the stellar radius that results from the measured sink ac- 
cretion rates using the stellar radius model shown in Figure [T] The 
expected trend of an increasing radius which reaches a sharp peak 
and then rapidly decreases is still found. However, variation due to 
the dependency upon accretion rate is now superimposed on top of 
this trend. Our model is semi-analytic, not a real stellar evolution 
model, and therefore the sharpness in the variation is most likely 
artificial. In reality, the protostar would only be able to respond 
to changes in the accretion rate according to its Kelvin-Helmholtz 
time. However, allowing the radius to vary along with the accretion 
rate decreases the variation in the luminosity, which is proportional 
to M/R t , and hence this represents a conservative choice for our 
purposes. Figure [10] also shows the accretion luminosity from the 
sink over the same period. Initially the variation in the accretion 
rate dominates both the stellar radius and the luminosity. However, 
over time the variation in the stellar radius becomes smaller as the 
star leaves the adiabatic phase and starts steadily contracting. Dur- 
ing this later phase the actual mass of the star is no longer changing 
so rapidly and there is a clear trend in the radius. Consequently as 
L a cc = GM^M /R t , the luminosity is no longer so noisy. 



5 DISCUSSION 

The effect of accretion luminosity feedback is a modifying factor 
affecting fragmentation in minihalos, rather than a dominant one. 
Inter-halo variability produces a greater variation in the number of 
fragments formed than feedback effects, and the number of frag- 
ments chiefly depended on the initial conditions of the halo which 
was being considered. In halos in which a large number of frag- 
ments formed, N-body interactions and ejections are at least as im- 
portant as accretion luminosity and play an important role in the ac- 
cretion histories of the protostars. It was only for more slowly frag- 
menting halos in which only a few fragments were initially formed, 
that the accretion feedback became important. At this point, the 
chief role of the accretion feedback was to delay fragmentation. 
This limits the number of fragments that can form before the largest 
sink particle becomes massive enough to start emitting significant 
quantities of ionising radiation. Even with the inclusion of feed- 
back, a small cluster of sink particles was formed in all the miniha- 
los studied here. 

However the analysis of the accretion rates and dynamics has 
raised some interesting questions in its own right. The accretion 
rates onto the sink particles were highly variable, even periodically 
so in some cases. The accretion rates were recorded at the sink 
accretion radius of 20 AU instead of directly onto a protostar, and so 
some of this variation may be smoothed by an inner disk. However, 
it is unlikely this effe ct could be removed e ntirely. Population III 
stellar modeling (e.g. lOmuk ai & Palla 2003) typically assumes a 
constant accretion rate, and it is unclear how a variable accretion 
rate would affect the stellar evoluti on of Popul a tion II I protostars. 

Moreover, as highlighted by lGreifetal.l ( l20T 1), there are a 



large number of encounters that can produce ejections. In treating 
encounters we use a simple sink particle prescription that consid- 
ers the stars as point masses, whereas in reality the protostars are 
extended gaseous objects with radii of typically around 100 Rm. 
Tidal forces will be strong during such an encounter raising the 
possibility the protostars might merge. The two approaches that are 
typically used for close encounters with sink particles are either to 
model the interaction as if it were occurring between point masses 
(as done here) or to combine the particles together as 'sticky' par- 
ticles. Both approaches are probably oversimplifications of the true 
picture, and given that we have found that encounters are more im- 
portant than feedback in determining the early history of population 
III protostars, it would be useful if this issue were addressed in the 
future. 

Another consequence of this work is the implications for pro- 
ducing 'dark stars'. It has been proposed that dark matter annihila- 
tions at the centre of a minihalo, where the dark matter density is 
at a maximum, could be a significant source of energ y that could 
support primordial stars with r adii of up to 10 AU (SpolyareTal] 
l2008ll2009l : lFreese et al.ll2008l) . However in the minihalos studied 
here, the protostars never precisely remain at the centre of their 
dark matter halos throughout this early stage of their evolution. It 
is, however, too early to say whether these findings exclude dark 
stars as the predicted dark matter annihilation luminosity is at least 
an order of magnitude greater tha n that found here from accretion 
feedback. iRipamonti et al.l found that when full gas chem- 

istry was include in a ID calculation, annihilation feedback was not 
sufficient to halt collapse and a normal population III protostar was 
formed. However, the contribution from annihilation feedback may 
be enough to prevent fragmentation, in which case the protostars 
would have no interactions and remain in the centre of their halos. 
This is a question that we are currently addressing. 



6 CONCLUSIONS 

We have introduced a prescription for heating from accretion 
luminosity into a re-simulation of five minihalos from cosmolog- 
ical initial conditions. We followed the evolution of these halos 
with and without feedback up until the point at which ionisation 
feedback would become significant. Our findings are: 

(i) Accretion luminosity delays fragmentation but cannot pre- 
vent it. 

(ii) The intrinsic variation in halo properties due to differences 
in their formation history generally has a larger effect on the num- 
ber of fragments formed than the accretion luminosity does. 

(iii) Halos in which a large number of fragments form rapidly 
are dominated by dynamical effects. It is only in more slowly frag- 
menting cases that form fewer fragments that accretion luminosity 
becomes effective. 

(iv) Accretion luminosity has little to no effect on the accretion 
rates of the protostars. On the other hand, dynamical ejections are 
an effective means of halting further accretion. 

(v) The accretion rates measured for the sink particles are highly 
variable and are quite different from the constant or slowly vary- 
ing accretion rates assumed in most pre-main sequence models for 
Population III stars. 
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Figure 9. The evolution of the sink particles in halos 1-4 with feedback. The most massive sink is shown in red. The accretion rates onto the sink particles are 
highly variable due to N-body dynamics causing sinks to orbit through their accreting gas streams. 
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Figure 10. Estimated stellar radius (solid line) and luminosity (dotted line) 
during the evolution of a typical massive sink. 
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